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Exact analytical expressions for the matrix elements of the Uehling potential in a basis of explic- 
itly correlated exponential wave functions are presented. The obtained formulas are then used to 
compute with an improved accuracy the vacuum polarization correction to the binding energy of 
muonic and pionic molecules, both in a first-order perturbative treatment and in a nonperturbative 
approach. The first resonant states lying below the n — 2 threshold are also studied, by means of 
the stabilization method with a real dilatation parameter. 
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The involvement of muonic molecular ions in nuclear fusion as fusion catalysts, through the Vesman mechanism 
generated great interest for precise energy level calculations in small muonic molecules 0. In particular, precise 
i— |, knowledge of the binding energy of the weakly bound state (L = l,v = 1) in ddfi and dtfi is required to predict the 
4m ' temperature dependance of molecular formation rates. The analysis of ddfi fusion experiments performed at PNPI 
actually resulted in a very precise determination of the (L = 1, v = 1) binding energy (with 0.7 meV uncertainty), in 
impressive agreement with theory [3| . Knowledge of the spectrum of resonant states below the n — 2 threshold is also 
useful for evaluating their impact in the muon catalyzed fusion cycle [3, Q • 

Exotic molecular ions also play a role in the interpretation of spectroscopy experiments in muonic or pionic atoms. 
The existence of \xp atoms in the metastable 2S state, a prerequisite for the measurement of the 2S-2P Lamb shift @ 
was observed through a quenching mechanism by collisions with H2 which involves resonant states of ppfi below the 
n = 2 threshold (tJ. In experiments on pionic hydrogen or deuterium Q, atoms are produced from highly excited 
states through an atomic cascade in which resonances of ppw or ddn may be populated @; the properties of these 
resonances are useful input parameters for an accurate modeling of the atomic cascade, which is indispensable to 
understand the observed line shape and extract the strong interaction broadening. 

Some of these applications (most notably muon catalyzed fusion studies) require accurate energy level calculations, 
which means that leading corrections to the nonrelativistic energies have to be taken into account. In muonic systems, 
by far the largest correction originates from the vacuum polarization contribution to the interaction energy, whereas 
in pionic systems the strong interaction shift is of the same order S. The first-order polarization correction to the 
interaction potential is usually referred to as the Uehling potential it is given by a nonelementary integral over 
a parameter. Most calculations of the Uehling correction in three-body systems have been performed by means of a 
numerical integration of its matrix elements, either with a Gaussian [5|, [3, [H| or exponential (l2l . [l3j basis set. An 
analytical expression of its matrix elements in a correlated exponential basis set was published in Ref. However, 
that expression is quite complicated, and numerical results obtained from it [l5[ are in disagreement with those of 
other authors. 

In the present work, we give in Sec. |TI]a more compact analytic expression for the matrix elements of the Uehling 
potential in a correlated exponential basis set, which greatly simplifies its application in actual calculations. These 
results may also be applied to calculations with the generalized Hylleraas expansion 16;]. The obtained expressions 
are then used in Sec. IHII to obtain a new set of reference results for bound and resonant state energies in muonic and 
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II. MATRIX ELEMENTS OF THE UEHLING POTENTIAL 



We use atomic units, scaled to the mass m of the lightest particle of the studied three-body system (e.g. the muon 
mass in the case of muonic molecules) . The Uehling potential between two particles of charges Z\ , Z 2 reads 



Mr) = f due— ^E1!&±A (!) 
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with x — (at sc m) 1 (here at sc represents the fine structure constant). We consider a variational expansion of the 
three-body wavefunction in the form 
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*(rx,r 2 ,r 12 ) = £ C„ e — -^-/J-ni^-n, ^(r^r,), (2) 

71=1 

where n, r 2 , n 2 are the interparticle distances and (r 1, f 2 ) are bipolar spherical harmonics. a n , /3 n , j n are real 
exponents satisfying the relations a n + j3 n > 0, a n + j n > 0, and (3 n + j n > 0. The matrix elements of the Uehling 
potential in such a basis set involve integrals of the form 



m fff j r f 00 Vu 2 - 1 (2u 2 + l) 
A 7) = / / d ri dr 2 dr 12 r[rTr^e- ar ^^-^ / aW 2 ""'*- ^ i 



(3) 

where r« = ri, r 2 , and ri 2 for T4p K>p (^2) and V vp {r\ 2 ) respectively, and l,m,n are non- negative integers. These 
integrals can be generated from Zo,o,o( Q! ! P> 7) by partial differentiation with respect to a, /?, and 7, as it is usually 
done in the case of the Coulomb potential (see e.g. [13, EH). The basic integral to be calculated is thus 



J$,o(a»&7) = /// dndr 2 dr 12 e -«i-**-T-« J°° due -^n ^^ 1 \* u + 1 ) (4) 

The first step is to change the order in which the integrations over space coordinates and the parameter u are 
performed. For V vp (ri) one obtains 



The integral over spatial coordinates is well known [17|, [3 an d reads 2/(/3 + 7) (a + /3 + 2xu)(a + 7 + 2a;w), so that 

I§ fi (a,M)= 2{ pl i)x2 h{a,b), (6) 
where a = (a + /3)/2x, b = (a + j)/2x, and 



, /- 00 V« 2 - 1 (2u 2 + 1) 
h(a,b)= du——r: \- -i. 7 

V ' 7 A w 4 (u + a)(u + 6) 1 ; 

The integral ([7]) can be obtained analytically by standard procedures (the work can be done using a symbolic com- 
putation program such as Mathematica): 

37r(a + b) [2 (a 2 + b 2 ) + 3a 2 b 2 ] - ab [12 (a 2 + ab + b 2 ) + 20a 2 6 2 ] 
= 



\/i - a 2 (l + 2a 2 ) arccos(a) \/I - fe 2 (l + 2o 2 ) arccos(o) 
+ a 4 (a-o) b 4 (a~b) ' (8) 

Since the last two terms in this expression diverge for a — b, one should study the limit b — > a. The result is 

, N 3tt (4 + 3a 2 ) -2a (12 + 11a 2 ) (2a 4 - a 2 - 4) arccos(a) 
h(a, a) = — ^ >— -g-i ^ + A r7== — ■ 9 

6a 5 a 5 Vl - a 2 

For 5* states, the matrix elements of V vp (ri) involve the integral 

r (i) , fl , 5 2 4S, (tt,^7) nn . 

4ii(«.A'r) = ^ (10) 
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Straightforward (but tedious) algebraic manipulations lead to 

h(a,b) h(a,b) h{a,b) 



(P + j)x 2 



(11) 



where 

h{a,b) 



(/3 + 7 ) 2 4x(fi + 7) 8a; 2 
3tt [4 (a 4 + a 3 b + a 2 b 2 + ab 3 + b 4 ) + 3 (a 4 6 2 + a 3 b 3 + a 2 & 4 )] - 2ab{a + b) [12 (a 2 + b 2 ) + lla 2 b 2 ] 

(12 



6a 5 o 5 

(4 + a 2 - 2a 4 ) arccos(a) (4 + b 2 - 2o 4 ) arccos(b) 



h(a,b) 



a 5 {a-b)Vl-a 2 b 5 (a - fyy/l^P 

37<a + b) [4 (a 2 + o 2 ) + 2ab + 3aW\ - 2a& [l 2 (a 4 + o 4 ) + 1 1 (a V + a 2 ?* 4 ) - 6 (a 3 &+ a 2 o 2 + C& 3 ) - lOa 3 ^ 3 ] /{a - b) 2 

6a 5 b 5 

(46 - 6a + a 2 fe - 3a 3 - 2a 4 6 + 6a 5 ) arccos(a) (4a - 6b + ab 2 - 3o 3 - 2ao 4 + 6b 5 ) arccos(fr) 
+ a 5 {a-b) 3 VT^ b^a-bfVT^W ' (13) 



In the case a = b, these expressions are replaced by 

3tt (20 - 11a 2 - 9a 4 ) - 2a (60 - 23a 2 - 28a 4 ) (20 - 21a 2 - 6a 4 + 4a 6 ) arccos(a) 

'*"> = - JiT^> U ("■> 

3 7 r(20+6a 2 )(l-a 2 ) 2 -a(l20-184a 2 + 23a 4 + 32a 6 ) (40-88a 2 +45a 4 + 10a 6 -4a 8 )arccos(a) 

I^ia.a) = ^ ^-rz . (14b) 

6a? (1- a 2 ) 2 2a7(l-a 2 ) 5/2 

The matrix elements of V vp (r 2 ) (resp. V vp (ri2)) can be deduced from this result by interchange of the parameters a 
and j3 (resp. a and 7). 

For P states, three integrals are needed: 1^1 1, IqI 1, and 3 . Their expressions are too lengthy to be reported 
here, but they can be easily evaluated by symbolic calculations programs and translated into C or FORTRAN code. 
For higher values of the orbital angular momentum, it is doubtful whether evaluation of analytical formulas remains 
advantageous with respect to numerical integration, because of growing calculation time and numerical instabilities. 

III. NUMERICAL APPROACH AND RESULTS 
A. Numerical approach 

In this Section, we present the results of variational calculations using the non relativistic three-body Hamiltonian 

H = --^-Vl i - -^V 2 2 - lv ri V r2 + (15) 

zmi 27712 m ri f2 ri2 

Here, the nuclei are numbered by 1 and 2, and the light particle (muon or pion) by 3. The notations n = 7*13, r 2 = 7'23 
are used. m\ and 777,2 are respectively the 1 — 3 and 2 — 3 reduced masses. The vacuum polarization correction to the 
binding energy is determined from first-order perturbation theory: 



AE™ = AE$ - ((V vp (n)) + (V vp (r 2 )) + (V vp (r 12 ))) (16) 



where AE$ is the first-order shift of the related atomic threshold. 

The Uehling potential behaves like ln(r)/r at r — > [19]. This not too singular behavior enables good convergence 
of a nonperturbative calculation, where the vacuum polarization potential V vp (ri) + V vp (r 2 ) + V vp (ri 2 ) is directly 
added to the Coulomb Hamiltonian H before diagonalization. The correction to the binding energy is then 

AE b = E% U) - E^ - [E$p - E^) (17) 

where E^ cu ) and E^ ){CU) are the energies of the three-body state and of the related atomic threshold, obtained 
with the Coulomb (C) or Coulomb+Uehling (CU) potential. While corrections beyond the first order are not useful 
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in themselves at the present level of theoretical accuracy, this provides a simple and reliable way of evaluating higher- 
order corrections and thus controlling the accuracy of the results. In addition, the perturbative approach fails for 
weakly bound resonant states close to a n > 2 atomic threshold. One way to understand this is to consider that the 
lifting of the atomic manifold degeneracy induced by the Uehling potential modifies the long range behavior of the 
atom- nucleus interaction potential, from a 1 A" 2 dipole potential to a 1 /r 4 induced dipole potential. A nonperturbative 
calculation is thus mandatory in such cases [t| [2l|] . In all the tables below, we give both the first-order perturbation 

result A£^ and higher-order corrections A£^ >1 ' ) = AE^ — AE^K 

The expansion ([2|) was used, with real exponents a n , (3 n and 7„ generated in a pseudorandom way in intervals 
[A\, A 2 ], [Bi,B 2 ] and [Ci,Cy respectively [22|, [Hj]. Here the variational parameters are the bounds of the intervals, 
and were optimized separately for each calculated level. Basis sets of N = 1000 — 2500 vectors were used to obtain 
good convergence of the results. 

It should be noted that complex exponents are generally better suited for molecular systems (24|. The analytical 
formulas of Sec. [Til are still valid for complex exponents a n , /3 n , j n , provided their real parts satisfy the relationships 
Re[a n + fi n ] > 0, Re[a„ + j n ] > 0, and Re[/3 n + j n ] > 0. However, with an expansion that uses complex exponents 
and/or using complex coordinate rotation [25j to study resonant states, numerical problems appear when the Uehling 
potential is included in the Hamiltonian. This suggests that the Uehling potential may not be dilation analytic 26] . 
A rigorous analysis of this point is beyond the scope of the present paper, but would certainly be useful for further 
studies with the Uehling potential. 

For nonperturbative calculations, it is important to add higher exponents in the basis set in order to describe 
accurately the behavior of the Uehling potential at small r. This is done by adding several subsets defined by 

- a AO) _ AO) 

. r n AO) An) _ n A0) ^> 

Typically r ~ 3 — 5, and n max = 1 — 2. We add similar basis sets for r 2 (if the basis is not symmetrized) and r\%. 

With the above-mentioned typical basis size, quadruple-precision arithmetic is generally required to maintain suffi- 
cient numerical stability. However, the derived expressions of the Uehling potential's matrix elements are numerically 
unstable (for a « 6), so that sextuple-precision arithmetic had to be used in most cases. For the weakly bound 
[L = 1, v = 1) states in ddfi and di/x, which require the largest basis sets, octuple precision proved necessary. 

We used the latest CODATA (2010) values [27| of the particle masses (muon, proton, deuteron and triton) and of 
the fine structure constant. For the pion mass, the latest value from the Particle Data Group [28|, was used. The 
quantity x appearing in the expression ([1]) of the Uehling potential is — 0.6627515411 for muonic systems, and 
Xn = 0.5017207015 for pionic systems. 




B. Results 



We first determined the vacuum polarization shift of the IS and 2S atomic thresholds, both in the perturbative 
and nonperturbative approaches, using a variational approach similar to the one described above. The radial atomic 
wavefunction ^f(r) is expanded on a set of N = 50 — 100 exponentials e~ Q '* r with pseudorandomly chosen real 
exponents. Results are summarized in Tabic [TJ 



Atom 


state 




AE 


MP 


IS 


-1.898 829 6 


-1.900 865 8 




25 


-0.219 584 


-0.219 737 2 


fid 


15 


-2.129 272 6 


-2.131642 2 




25 


-0.245 319 4 


-0.245 494 5 


fit 


15 


-2.214430 5 


-2.216 9261 




25 


-0.254 804 


-0.254 987 2 


Tip 


15 


-3.240 8019 


-3.244 916 5 




25 


-0.368 276 3 


-0.368 560 


nd 


15 


-3.732175 


-3.737120 2 




25 


-0.422 196 4 


-0.422 529 5 



TABLE I: Vacuum polarization shift of the 15 and 25 atomic states of muonic and pionic atoms, in eV. Both the first-order 
perturbation result AE*- 1 ' and the nonperturbative result AE are given. 
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Table [TTI gives the energies of all the bound states of muonic molecules with orbital angular momentum L = 0,1. 
The results are in perfect agreement with earlier calculations [29], with an accuracy improved from 0.1 meV to 1 
fj, eV. The contribution from higher perturbation orders is also obtained (for the first time to our knowledge), and 
typically amounts to a fraction of meV for the ground vibrational state. Precise experimental results are available 
only for the (L = l,v = 1) state of ddfi where there is good agreement with theoretical predictions [30Tl32l | that 
also take leading relativistic and nuclear structure corrections, as well as corrections caused by the finite size of the 
(dd/i)dee molecular complex. The discrepancy is only 0.5 meV, while experimental and theoretical uncertainties are 
respectively of 0.7 and 0.4 meV. The 0.097 meV difference (-8.657 meV instead of -8.56 meV) between our new result 
for the Uehling correction, and the value of does not alter the agreement with experimental data. The newly 
obtained contribution from higher perturbation orders (0.003 meV) is currently not relevant in view of the overall 
theoretical uncertainty. 



Molecule 


L 


V 




E b 














(cV) 


(moV) 


(meV) 


ppfi 








253 


.150 104 


284.875 


0.430 




1 





107 


.265 303 


50.581 


0.089 


pdn 








221 


.547 587 


234.419 


0.376 




1 





97 


.497 678 


21.445 


0.053 


ptn 








213 


.838 459 


222.385 


0.365 




1 





99 


.126 024 


21.009 


0.055 


ddfi 








325 


.070 580 


412.131 


0.657 







1 


32 


.844 224 


39.129 


0.074 




1 





226 


.679 812 


226.216 


0.358 




1 


1 


1 


.974 980 


-8.657 


0.003 


dt/i 








319 


.136 858 


402.275 


0.653 







1 


34 


.834 420 


28.074 


0.061 




1 





232 


.469 701 


233.597 


0.377 




1 


1 





.660 329 


-16.604 


0.013 


ttfj, 








362 


.906 436 


480.211 


0.781 







1 


83 


.770 686 


99.858 


0.172 




1 





107 


.265 303 


331.988 


0.534 




1 


1 


45. 


.205 712 


34.072 


0.072 



TABLE II: Vacuum polarization correction to the binding energies for bound states of muonic molecules, obtained using the 
variational expansion ([2]| with real exponents. The binding energy Eb calculated with the pure Coulomb potential is given in 
the fourth column. The vacuum polarization shift at first order of perturbation theory is given in the next column. The last 
column shows the difference between results of non-perturbative and first-order perturbative treatments. 

Results for the bound states of pionic molecules are given in Table Mil We have limited our study to ddir and 
ppir, which could play a role in the interpretation of pionic hydrogen and deuterium spectroscopy experiments [8J|. It 
should be noted that accuracy is much less essential than for muonic systems, because (i) experimental resolution is 
limited to about 10 /ieV by the pion lifetime r = 26 ns, and (ii) theoretical accuracy is limited to a fraction of meV 
by the 2.5 ppm relative uncertainty on the pion mass. However, the vacuum polarization correction is relevant since 
it typically amounts to a fraction of eV for the ground vibrational state. 



Molecule 


L 


V 


E b 












(cV) 


(mcV) 


(moV) 


ppir 








294.859 450 


431.020 


0.763 




1 





80.227 512 


6.808 


0.055 


dd% 








392.301211 


660.791 


1.237 







1 


15.777113 


19.426 


0.053 




1 





237.301 428 


291.614 


0.556 



TABLE III: Same as Table (TT] for bound states of the pionic molecules ppn and ddir. 



In the following, we consider quasibound states (or resonances). In view of the problems with complex coordinate 
rotation mentioned in Sec. MI Al we used the stabilization technique with a real dilatation parameter, similarly to 
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Ref. [20|. The accuracy of this method is limited by the width of the resonances. In the following tables, we report 
the widths calculated in Ref. [33[ in order to explain the accuracy of the results. While a complete investigation of the 
resonance spectrum would lie beyond the scope of this paper, we give illustrative results for the first two vibrational 
and rotational states below the 2S threshold. 

Among the muonic molecules, we have considered dd[i and dt/i, in which fusion research has been the most active 
(see Table IIV|) . The involvement of resonances was originally proposed in the framework of d — t fusion, whereas its 
impact in d — d fusion is expected to be much less important 0, [34|. In the case of dt\x, our results are in good 
agreement with those of Ref. [20l |. and represent an improvement in accuracy by 2-3 orders of magnitude. 



Molecule 


L 


V 


E b 


r 


AE? 


AE b 


AEb 








(oV) 


m 


(moV) 


(mcV) 
This work 


(mcV) 

m 


ddfi 








218.11160 


1.9 - 


-54.77 


-54.79 









1 


135.279 02 


5.8 - 


-82.76 


-82.79 






1 





211.924 50 


5.8 - 


-58.44 


-58.46 






1 


1 


130.350 1 


15.3 




-85.5 




dtfi 








217.889 86 


3.0 - 


-59.83 


-59.86 


-60 







1 


139.73140 


7.2 - 


-86.66 


-86.70 


-85 




1 





212.545 744 


0.5 - 


-63.006 


-63.030 


-63 




1 


1 


135.379 516 


0.9 - 


-89.069 


-89.104 


-91 



TABLE IV: Vacuum polarization correction to the binding energy for resonant states of the muonic molecules ddfi and dtfi 
below the n = 2 threshold, obtained using the variational expansion ([2]) with real exponents. The binding energy Eb obtained 
with the pure Coulomb potential is given in the fourth column. The fifth column contains the resonance widths taken from [33j], 
which give a measure of the precision of the results. The vacuum polarization shift, both at first order of perturbation theory 
and in a nonperturbative treatment, are given in the next two columns. The first-order result is given only in the cases where 
the precision is sufficient to evidence the difference with the nonperturbative result. 

Table |V] summarizes results for the pionic molecules ppir and ddir, where resonant states play a role in the deexci- 
tation cascade of pionic atoms @. In the case of ppir, our results are in good agreement with those of Ref. [9|, and 
bring an improvement in accuracy by 1-2 orders of magnitude. Note that the binding energies of the resonances we 
have studied are large enough for the perturbative approach to yield precise results. The difference with the result of 
a nonperturbative calculation is typically of 20-40 /ieV only. 



Molecule 


L 


V 


E b 


r 


AE^ 


AEb 


AE b 








(cV) 


( M eV) 

en 


(mcV) 


(mcV) 
This work 


(moV) 
1 


ppir 








236.173 


1.5 




-78 


-80 







1 


100.146 


1.9 




-136 


140 




1 





220.381 8 


0.20 




-92.0 


-90 




1 


1 


89.641 


0.38 




-145 


150 


ddir 








275.280 3 


0.050 




-85.5 









1 


156.8218 


0.097 




-139.7 






1 





265.180 84 


0.004 1 


-93.87 


-93.90 






1 


1 


149.088 74 


0.005 4 


-145.56 


-145.60 





TABLE V: Same as Table HVl for resonant states of the pionic molecules ppir and ddir below the n = 2 threshold. 

In conclusion, we have shown that the matrix elements of the Uehling potential in a basis of correlated exponential 
functions may be expressed in an analytical form. We have used the obtained expressions to calculate the vacuum 
polarization shift for a wide range of bound and resonant states in muonic and pionic molecules, either for the first 
time, or with a greatly improved accuracy. The excellent agreement with earlier calculations which used a numerical 
evaluation of matrix elements fully confirms the validity of the analytical formula. 
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